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ABSTRACT 


A  method  is  developed  in  terms  of  the  Cauer  Second  Form 
representation  of  continued  fractions  as  a  means  of  designing 
linear  single-input  single-output  (SISO)  control  systems. 
Optimal  closed  loop  solutions  corresponding  to  a  set  of 
prescribed  eigenvalues  are  obtained  through  minimization  of 
a  quadratic  performance  index.  The  partitioning  method  of 
the  Cauer  Second  Form  for  system  simplification  is  presented 
with  a  simplified  inversion  technique  for  the  reduced  order 
system. 
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I.  INTRODUCTION 


The  purpose  of  this  research  was  to  develop  an  algorithm 
for  obtaining  optimal  closed  loop  solutions  corresponding 
to  a  set  of  prescribed  eigenvalues  for  single-input  single¬ 
output  (SISO)  control  systems.  It  was  desired  that  the 
algorithm  be  adaptable  to  digital  computer  techniques  and 
unrestricted  by  system  order. 

The  Cauer  Second  Form  for  system  dynamics  representation 
was  chosen  over  other  alternatives  because  of  the  regular 
pattern  of  the  state  and  output  matrices,  and  the  method 
of  linear  system  simplification. 

In  Chapter  II,  several  basic  properties  of  both  Cauer 
First  and  Second  Forms  are  presented  from  the  theory  of 
continued  fractions.  A  simple  and  efficient  algorithm  is 
also  developed  for  inversion  of  the  continued  fraction  in 
either  form,  independent  of  Routh's  algorithm. 

In  Chapter  III,  the  method  of  linear  system  order 
reduction  based  on  the  Cauer  Second  Form  is  amplified. 

The  emphasis  on  this  area  was  primarily  to  elucidate  the 
various  methods  previously  employed  for  system  simplifi¬ 
cation  . 

The  original  theoretical  work  of  this  thesis  is  pre¬ 
sented  in  Chapter  IV.  The  objective  was  to  obtain  closed 
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loop  solutions  corresponding  to  a  prescribed  set  of  eigen¬ 
values.  While  minimizing  a  certain  cost  function,  which 
met  desired  system  characteristics.  It  is  shown,  by  examples 
that  the  derived  algorithm  is  equally  capable  of  handling 
systems  with  multiple  and/or  complex,  as  well  as,  unique 
sets  of  real  eigenvalues. 

The  final  chapter,  Chapter  V,  presents  a  discussion  of 
results  and  suggests  areas  for  future  study. 


PROPERTIES  OF  CAUER  FIRST  AND  SECOND  FORMS 


A.  CLOSED  LOOP  SYSTEM  IN  CAUER  FIRST  AND  CAUER  SECOND 
FORMS 


Consider  the  closed  loop  transfer  function  given  by: 


Y(s) 

uTST 


n-l 

I 

L=0 


n 


s 


+ 


b.  s1 
1 

n-l 

Z 

i  =  0 


a .  s 
1 


i 


(2-1 


with  block  diagram  as  given  in  Figure  2.1.  Equation  (2-1) 
can  be  expanded  into  the  Cauer  Forms  of  continued  fractions 
as  follows. 

1 .  Cauer  First  Form 

a.  Arrange  the  numerator  and  denominator  poly¬ 
nomials  in  descending  order. 


b. 


Perform  continued  division 


Programming  (Bush  Form) 


2.  Cauer  Second  Form 


a.  Invert  the  numerator  and  denominator  and  arrange 
the  polynomials  in  ascending  order. 


Y(s) 

uTiT 


a  +  a, s  +  _ a  0sn"2  +  a  ,sn_1  +  sn 

o _ 1 _ n- 1 _ n-1 _ 

bo  *  V  *  . tbn-2sn"2  *  Vl5"'1 


(2-4) 


b.  Perform  continued  division. 


Y(s) 

uTsT 


1 


1 


h3  + 


(2-5) 


or 


1 


hl  + 


h2  + 


h3  + 


h4  + 


(  2-6) 


Block  diagrams  of  both  systems  are  shown  in  Figures  2.2  and 


2.3. 
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Cauer  First  Form 
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B.  PHYSICAL  INTERPRETATION  OF  DOMINANT  TERMS 
RESULTING  FROM  CONTINUED  FRACTION  EXPANSION 

It  is  known  that  the  most  dominant  terms  in  equations 
(2-3)  and  (2-5)  are  the  first  quotients,  h^s  and  h^,  re¬ 
spectively.  A  meaningful  interpretation  for  these  terms 
can  be  accomplished  by  applying  the  initial  value  and  final 
value  theorems. +  Letting  Y(s)/U(s)  =  F(s),  by  an  asymptotic 
expansion  approximation: 


1.  For 

Cauer  First  Form 

lim 

f(t)  =  lim  sF( s )  = 

hl 

(2-7) 

t-*o 

S-H» 

and 


lim  f(t)  -  lim  sF(s)  ~  +  ^4  ‘ 

(2-8) 

t-*-°°  s-*-o 

2. 

For 

Cauer  ! 

Second 

Form 

lim 

f(t) 

~  lim 

sF(  s ) 

=  h2  +  h4 

(  2- 

9) 

t-*-o 

S-*®3 

lim 

f(t) 

z  lim 

sF(s) 

=  hr 

(  2- 

10) 

t-«o 

s-+o 

lim  f(t)  must  exist. 

15 


Equations  (2-7)  and  (2-10)  are  of  considerable  interest  since 
they  involve  the  dominant  term,  h^.  The  implication  is  that 
the  Cauer  First  Form  emphasizes  the  initial  or  transient 
response  of  the  system;  whereas,  the  Cauer  Second  Form  em¬ 
phasizes  the  final  or  steady  state  response  of  the  system. 

In  general,  the  quotients  lower  in  position  ip  the  continued 
fraction  expansion  have  less  influence  on  the  performance  of 
the  system  as  a  whole,  (h ^  has  less  influence  than  h^,  where 
i<j).  Because  many  systems  must  meet  a  set  of  steady  state 
conditions,  the  Cauer  Second  Form  will  be  used  for  the 
prescribed  eigenvalue  problem. 

C.  CONTINUED  FRACTION  INVERSION 

The  theory  of  continued  fractions  was  first  associated 
with  Routh's  Algorithm  by  Wall  in  1945,  [1]  and  [2].  The 
following  year  Frank  [3]  extended  and  modified  Wall's  work 
to  include  complex  coefficients.  Both,  however,  applied 
Routh's  algorithm  only  to  continued  fraction  expansions, 
not  to  the  problem  of  inversion. 

In  1969,  Chen  and  Shieh  [4]  developed  an  algorithm 
method  for  converting  a  continued  fraction  into  a  rational 
fraction  of  two  polynomials.  Their  method,  which  makes  use 
of  Routh's  algorithm,  is  presented  below. 

If  the  elements,  h^,  are  known  for  any  continued 
fraction,  then  the  state  and  output  equations  can  be 
written  immediately  from  Figures  2.2  or  2.3. 
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The  coefficients  of  the  characteristic  polynomial, 

|sI-H],  become  the  first  row  elements  of  the  required  Routh 
array.  The  next  sequence  of  steps  in  determining  the  (j+2)th 
row  of  the  Routh  array  is  to  successively  let 


h.  =  0 
1 

and 

h.+1  =  0  (2-15) 

for  j s[l ,  3 , 5 , . . . 2n_l] 

and  evaluating  the  remaining  (n-k)x(n-k)  ,rJi"  matrix,  where 
k  =  (j+l)/2,  up  to  k=n-l,  i.e.,  for  n  arbitrary  and  j=l; 
the  3rd  row  of  the  Routh  array  becomes  (after  h^  and  h3  are 
set  equal  to  zero)  the  coefficients  of 


sI-Hfl , 


where  the  (n-l)x(n-l)  matrix  is: 


h4h3  h6h3 


>h2nh3 


h4h3  h6(h3+h5)  . h2n(h3+h55 


h4h3  hg(h3+hj) 


,h2n(h3+"  *h2n-l) 


(2-16) 


This  process  repeats  until  the  system  state  matrix  is 
reduced  to  a  single  element,  ]_»  yielding  the  (2n-l)th 

row  in  Routh 's  array.  It  is  observed  that  each  successive 
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odd  numbered  row  contains  one  less  element  than  it's 


predecessor.  By  inserting  leading  zeros  in  the  3rd,  5th, 
(2n+l)th  row,  the  matrix,  P,  is  formed. 


3rd 

Pll 

P12 

P13  . 

. 1 

5th 

0 

P22 

P23  . 

. 1 

7th 

0 

0 

P33  . 

. 1 

(2-17) 


(  2n+l  >thl 


0 


0  0 


1 


The  matrix,  P,  is  the  linear  transformation  matrix  required 
to  obtain  a  linear  system  in  Cauer  Second  Form  from  phase 
variable  (canonical)  form.  Continuing,  the  second  row  of 
the  Routh  array  is  obtained  from  the  output  matrix,  L, 
and  the  above  transformation: 


c  =  L.z  (continued  fractions) 

z  =  Px  (linear  transformation) 

y  =  Cx  (output  equation,  phase  variable  form) 

Therefore , 

C  =  L  P  .  (2-18) 

C  is  an  (lxn)  vector  whose  elements  are  the  second  row  of 
the  Routh  array. 

Consider  the  Routh  array  as  an  (n+l)x(2n+l)  matrix  with 
typical  element  .  The  quotients,  hu  ,  of  the  continued 
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fraction  expansion  can  be  expressed  as: 


il 


ri+l ,1 


(2-19) 


From  this  relationship  and  knowledge  of  how  the  Routh 
array  is  generated,  the  remaining  even  numbered  rows  of  the 
array  can  be  found.  The  transfer  function  as  a  ratio  of  two 
polynomials  is  written  as: 


T(s) 


n 

Z 

ill 

n+l 

Z 

L  =  1 


(2-20) 


Chen  and  Shieh  [4]  contend  that  this  method  is  the 
easiest  in  attaining  the  inversion.  The  author  disagrees 
and  presents  a  simpler  iterative  method  based  on  the  in¬ 
version  technique  for  the  Generalized  Cauer  Form  given  by 
Goldman  [5].  The  method  is  equally  suited  to  both  Cauer 
First  and  Cauer  Second  Forms ,  requiring  no  prior  knowledge 
of  Routh' s  algorithm.  Assuming  all  h^'s  are  known,  and 
non-zero,  in  equation  (2-3)  or  (2-5),  let: 


ai  =  h2i-l 


(2-21) 


b. 


i 


(2-22) 


for  i  e[1 , 2 , . . . . ,n] . 
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1 .  Inversion  of  Cauer  First  Form 

Initialize  two  (n+lxl)  vectors  C  and  D. 


c  = 

Cc0  c1  c2  . . , 

**•  Cn] 

(2- 

23) 

D  = 

td0  di  d2  * '  ' 

...  dn: 

(2- 

24) 

dn-(i+l)+j  “  an-i  x  Cn-i+j  +  dn-(i+l)+j’ 


(2-29) 


where  j r[ 0 , 1 , 2 , . . . ,i]  are  substituted  in  ascending  order,  and 
(2-28)  is  solved  before  (2-29)  for  each  value  of  j.  Now, 
let  i=2  in  equations  (2-28)  and  (2-29)  and  repeat  the  same 
procedure.  The  index  "i"  is  incremented  until  i=n-l,  and 
(2-28)  and  (2-29)  are  solved  as  before  over  the  appropriate 
range  of  the  index  "j".  The  final  vectors,  C  and  D,  contain 
elements  which  are  the  coefficients  of  the  numerator  and 
denominator  polynomials,  respectively,  of  the  transfer 
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function  (or  driving  point  impedance  function): 


T(s) 


n 

E 

i  =  l 


C .  s 
x 


n-i 


n 


E 

j  =  0 


d. 

3 


n- 


s 


j 


Example : 


(2-30) 


_  10s2  +  171s+360 
T(s)  =  — - j - 

s  +71s  +702S+720 


(2-31) 


By  continued  fraction  division: 


10s2+171s+360 


P  +  71s2  +  702s  +  720^  .Is 
s3+17.1s2+36s 

53 . 8s2  +  666s+720 


53.9s2  +  666s  +  720  (l0s2  +  171s+36o(  1  .1855 

10s2+123.562s+133.58 
47. 438 s  + 226 .42 


47 .438S  +  226 .42  Is 3 . 9s  2+  66  6s  +  7 2 (A  1.1362s 

53.9s2+257. 258s 
408.742S+720 


408.742S+720  |47 . 438s+226 . .115 

47.438S+82.79 

143.63  (2-32) 
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* 


143.63 


r40  8 . 742s  +  720<.  2.8458s 

408 . 742s 

T2Q 


720  J14  3 .6  3(  .1995 

'143.63 
0 


Therefore,  the  transfer  function  in  the  form  of  equation 
(2-3)  is: 


T(s) 


,1s  + 


.1855  + 


1.1362s  + 


115  + 


2.8458s  +  1 


.1995 


(2-32) 


with 

hx  =  .1 
h2  =  .1855 
h3  =  1.1362 


h4  =  .115 

h5  =  2.8458 

h.  =  .1995.  (2-33) 

b 


How,  using  equations  (2-21)  and  (2-22); 


i 

i 
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MM i iMMtr 


.1855 


b2  =  h4  =  .115 
b3  =  hg  =  .1995. 


From  equations  (2-25)  through  (2-27), 


(2-34) 


c3  =  b3  =  .1995 

d2  =  a3xc3  =  2 .8458( .1995)  =  .568  (2-35) 

d3  =  1. 

Substituting  i=l  into  equations  (2-28)  and  (2-29)  for  j=0: 


c2  =  b2xd2+c2  =  .115(.568)  +  0  =  .0653 
dl  =  a2xc2+dl  =  1-1362( .0653)  +  0  =  .0742 
for  j  =  1: 


c  3  =  b2cd3  +  c3  =  .115(1)  +  .1995  =  .  3145 

d2  =  a2xc3+d2  =  1.1362( . 3145)  +  .568  =  .9353 

At  this  point,  j=i,  therefore  increment  index 
for  i=2,  j=0: 


c1  =  b1xd1+c1  =  .  185 5 (  .0742)  +  0  =  .0138 
dQ  =  a-Lxc1+dQ  =  .  1(  .  0138  )  +  0  =  .  00138 
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for  j=l: 

C2  =  ^]_x<^2+c2  =  *^55(.  9353)  +  .  0653  =  .  2388 
d1  =  a1xc2+d1  =  .  1 ( .  2 3 8 8  )  +  .  0742  =  .  0981 
for  j=2: 

c  3  =  b1xd3+c3  =  .  1855(  1)  +  .  3145  =  .  5000 

^2  =  alXC3  +  <^2  =  +  *9353  =  .  9853 

Now,  at  the  point  where  j=i,  and  i=n-l,  the  transfer 
function  is: 


m,_x  _  .0138s2+.2388s+.5  „ 

Tts)  -  - = - 2 -  (2-3 

.  0013  8s  +.  09 8  Is  +.  98  5 3s  +  l 

Multiplying  numerator  and  denominator  by  l/dQ  =  720  yields: 

T(s)  =  1°32*171B*360  ( 2-3 

s  +71s  +702S+720 


2 .  Inversion  of  Cauer  Second  Form 

Initialize  the  two  (n+lxl)  vectors  C  and  D,  (2-23) 
and  (2-24),  to  all  zeros,  except: 


= 

b 

( 2-38) 

n 

n 

=  1 

( 2-39) 

n-1 

_ 

a  xc 

( 2-40) 

n 

n  n 
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The  following  set  of  equations  are  first  solved  for  i=l. 


c  ...  =  b  .  x  d  ...  +  a  .  ,  .  . , 
n-l+j  n=i  n-i+j  n-i+j+1 


^n-(i+l)+j  ”  an-i  x  ^n-d+D+j  +  ^n-i+1 


(2-41) 

(2-42) 


where  jE [0 ,1 , 2 i)  are  substituted  in  ascending  order, 
and  (2-41)  is  solved  before  (2-42)  for  each  value  of  j. 


Next,  find  d^  according  to: 


d  =  a  .  x  c  , 
n  n-x  n 


( 2-43) 


Now,  let  i=2  in  equations  (2-41)  and  (2-42)  and  repeat  the 
same  recursive  procedure.  The  index  "i"  is  incremented  by 
one  until  i  =  n-1,  and  for  each  value  of  i,  (2-41)  and  (2-42) 
are  iteratively  solved  over  the  appropriate  range  of  the 
index  "j".  The  resulting  elements  of  C  and  D  are  the 
coefficients  of  the  numerator  and  denominator  polynomials 
of  the  transfer  function: 


T(s)  = 


n 

Z  e.  s 
L=1  1 


n-i 


n 


(2-30) 


Z  d .  s 
j  =  0  1 


n-1 


Using  the  sample  example,  (2-31); 


T(s) 
+for  j=i, 


10s2+171s+360 

s3+71s2+702s+720 

cn-i+ j  +1  =  0. 
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Place  the  numerator  and  denominator  terms  in  ascending 
order,  invert,  and  perform  continued  fraction  division: 


360+171s+10s' 


/720  +  702s+71s2  +  s3t 
720+342s+20s2 

360s+51s2+s3 


360s+51s+s" 


/360+171s  +  10s^ 
360+51s+s2 

120s+9s2 


1/s 


120s+9s ‘ 


f360s  +  51s  ~  +  s 

360s+27s2 

2  3 

24s  +sa 


24s  +  s  (120  +  9  st  5/s 

120+5s 


45 


4  ('24+ s(  6 

24 


S  (4i  4/s 

4 
0 


The  transfer  function  (2-31),  in  the  form  of  equation  (2 
is : 


T(s)  =  _ 1 


2 
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■Mintfitii 


2-44) 


-6) 


5  + 


s 


6  + 


s 


4 


(2-45) 


where 

\  =  2  \  -  S 

h2  =  1  h5  =  6 

h3  =  3  h6  =  4  (2-46) 

For  the  inversion  process,  using  equations  (2-21),  (2-22) 
and  ( 2-46) ; 

al  =  hi  =  2  bl  =  h2  =  1 

a2  =  h3  =  3  b2  =  h4  =  5 

a3  =  h5  a  6  b3  =  hg  =  4  (2-47) 

a.nd  from  equations  (2-38)  through  (2-40), 


G3 


=  4 


d3  =  a3  x  C3  =  6(4)  =  24.  (2-48) 

Now,  substituting  i  =  1  into  equation  (2-41)  and  (2-42),  for 
3  =  0: 


c2  =  b2  +  c3  =  5  +  4  =  9 


di  =  aj  x  ci  +  d2  =  3(0)  +1=1 
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(2-49) 


for  j  =1 ; 

03  =  b2  x  d3  +  0  =  5(24)  =  120 

d2  =  a2  X  t2  +  d3  =  3(95  +  24  =  51 

and  from  equation  (2-43): 

d3  =  a2  x  c3  =  3(120)  =  360 
for  i=2 ,  j=0: 

C1  =  bf  =  c2  =  1  +  9  =  10 

dQ  =  ax  x  cQ  +  d±  =  2(0)  +1=1 

for  j  =1  *, 

c2  =  b1  x  d2  +  c3  =  1(51)  +  120  =  171 

dl  =  al  x  ei  +  d2  =  2(10  +  51  =  71 

for  j  =2 ; 

d3  =  bx  x  d3  +  0  =  1(360)  +  0  =  360 

d2  =  a±  x  c2  +  d3  =  2(171)  +  360  =  702 

and  from  equation  (2-43): 

d3  =  a1  x  c3  =  2(360)  =  720. 
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Therefore , 


c  = 

Co 

10 

171 

360] 

it) 

II 

Cl 

71 

702 

720], 

(2-51) 

and  the  transfer  function  realized  from  equation  (2-30)  is: 


T(s) 


10s2+171s+36Q 

s3+71s2+702s+720 


(2-52) 


which  is  the  same  as  (2-31). 

This  completes  the  development  of  the  continued  fraction 
inversion  algorithms  from  Cauer  First  and  Second  Forms.  This 
iterative  procedure  is  easily  seen  to  be  computationally 
much  simpler  than  Chen  and  Shieh’s  method.  First,  it  does 
not  require  the  need  to  find  the  H  matrix,  (2-11);  and 
second,  it  does  not  necessitate  finding  the  coefficients 
of  n  characteristic  polynomials  of  diminishing  order.  This 
method  is  solely  based  on  equation  (2-6),  enumerating  the 
inversion  from  bottom  to  top.  As  by-product,  the  entire 
Routh  array  appears  in  the  intermediate  steps  as  can  be 
seen  from  the  Cauer  Second  Form  example: 


d3 

d2 

dl  d0 

720 

702 

71 

*3 

c2 

C1 

360 

171 

10 

d3 

d2 

dl 

360 

51 

1 

=3 

C2 

120 

9 

d. 

24 

1 
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where  rows  2(n-i)l  and  2(n-i)  are  taken  from  the  ith  iteration, 


i:e[0 ,1, 2 , .  . .  ,n-l]  .  "i  =  0"  implies  the  rows  come  from  the 
initialization  of  C  and  D.  The  last  row,  the  (2n+l)th,  is 
always  the  single  element  one. 

If  a  comparison  is  made  between  equations  (2-20)  and 
(2-30),  it  is  observed  that: 


c. 

i 


=  r2 ,n-i+l ’ 


ietl, 2 , . . . ,n] 


d •  ■  r^  y-j+2*  ] jl , 2 , . . • ,n3  ,  (2  —54 ) 

where  the  c^'s  and  d^'s  are  taken  from  the  (n-l)th  intera- 
tion  under  the  index  "i". 

It  is  also  observed  that  if  the  quotients,  h^'s, 
resulting  from  expansion  into  Cauer  Second  (First)  Form 
are  used  in  the  inversion  algorithm  presented  for  Cauer 
First  (Second)  Form,  then  the  c^'s  and  d^'s  in  the  (n-l)th 
iteration  represent  the  transfer  function  coefficients  in 
reverse  order.  This  is  shown  using  the  preceding  example. 
From  equation  (2-46); 


hx  =  2  h4  =  5 
h2  =  1  h5  =  6 
h3  =  3  h6  =  4 

and  equation  (2-47); 

a^  =  h^  =  2  b^  =  h2  =  1 
a  2  =  h  3  =  3  b2  =  h4  =  5 
a3  =  h5  =  6  b3  =  h6  =  4 
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Now,  using  the  inversion  scheme  for  Cauer  First  Form,  from 
(2-25)  through  (2-27); 

C3  =  b3  =  4 

d2  =  a3xc3  =  6(4)  =  24  (2-55) 

d3  = 

Making  the  substitution,  i=l,  in  equations  (2-28)  and  (2-29), 
for  j  =  0; 

c2  =  b2  x  d2  +  C2  =  5(24)  +  0  =  120 

=  a2  x  C2  +  d±  =  3(120)  +  0  =  360 

for  j=l,  ( j  =  i)  ; 

C3  =  b2  c  d3  +  d3  =  5(1)  +4=9 

d2  =  a2  x  c3  +  d2  =  3(9)  +  24  =  51 

for  i  =  2 ,  j  =  0; 

C1  =  bx  x  &  +  c1  =  1(360)  +  0  =  360 

dQ  =  aL  x  e  +  dQ  =  2(360)  +  0  =  720 

for  j  =1 ; 

c2  =  b1  x  d2  +  c2  =  1(51)  +  120  =  171 

d2  =  a2  x  c2  +  d2  =  2(171  +  360  =  702 
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for  j  =  2,  ( j  =  i) ; 


c3  =  b1  x  d3  +  c3  =  1(1)  +  9  =  10 

d2  =  a]_  x  c3  +  d2  =  2(10)  +  51  =  71 

and  dj  remains  unchanged,  equal  to  1. 
Therefore; 


•C  = 

[0 

360 

171 

10] 

D  = 

[720 

702 

71 

1], 

(2-56) 

and  the 

transfer 

function 

should 

be ; 

T(s) 


360s2+171s+10 
720s3+702s2+71s+1  . 


(  2-57  > 


Since  the  h^ ' s  from  the  Cauer  Second  Form  were  used  in  the 
inversion  algorithm  from  the  Cauer  Frist  Form,  the  vectors 
C  and  D  require  reversing  non- zero  elements,  resulting  in: 


[0 

10 

171 

360] 

[1 

71 

702 

720]  , 

(2-51) 

and  the  correct  transfer  function  is 

! 


_  10s2  +  171s  +  360 

US  J  -  — 5 - k — 1 - 

s  +71s  +702S+720 


(2-31) 
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A  digital  computer  program  (FORTRAN  IV)  has  been 
written  for  both  Cauer  First  and  Cauer  Second  Forms  and  is 
included  as  Appendix  3  with  documentation. 
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III.  LINEAR  SYSTEM  ORDER  REDUCTION  VIA  PARTITION 
OF  THE  CAUER  SECOND  FORM 


A  control  system, in  general,  can  consist  of  many  tens 
or  hundreds  of  elements.  In  such  cases,  the  problems  facing 
the  engineer  include:  (1)  too  many  variables  to  efficiently 
handle;  (2)  the  dimension  of  the  system  is  too  high  to  com¬ 
prehend;  and  (3)  the  modifications  needed  to  meet  required 
design  characteristics  are  difficult  to  ascertain.  A  logi¬ 
cal  approach  is  to  seek  procedural  methods  which  reduce  the 
order  of  the  system  to  a  manageable  size  yet  maintain  the 
basic  characteristics  of  the  full  dimension  model. 

A  number  of  different  methods  for  system  simplification 
have  been  proposed  for  the  reduction  of  high  order  dynamic 
systems  to  low  order  models  of  a  more  computationally  or 
analytically  tractable  nature.  The  approaches  used  are 
quite  different,  but  appear  to  fall  into  three  main  groups. 
The  first  is  to  ignore  those  modes  of  the  original  system 
which  contribute  little  to  the  overall  response.  Davison 
[6]  chose  to  neglect  eigenvalues  of  the  original  system 
which  are  farthest  from  the  origin,  retaining  only  the  do¬ 
minant  eigenvalues  and  hence  dominant  time  constants  in  the 
reduced  model.  The  shortcoming  of  this  technique  is  that 
many  systems  do  not  have  any  "dominant"  roots  [7]. 

Chidambara  [8]  essentially  finds  a  reduced  forcing  function 


35 


so  that  the  steady  state  values  of  the  lower  order  model 
agree  with  those  of  the  original  system.  The  consequence  of 
this  method  is  that  the  approximate  model  give  correct 
steady  state  values  but  incorrect  time  responses  because 
the  reduced  forcing  function  does  not  excite  the  modes  of 
the  two  systems  in  the  same  proportions  [6].  Marshall  [9] 
proposed  the  reduction  of  the  state  matrix  by  partitioning 
it  and  setting  certain  rate  variables  equal  to  zero  in  order 
to  maintain  the  original  steady  state  values.  This  techni¬ 
que,  like  Davison's,  is  based  on  dominant  roots  and,  there¬ 
fore,  exhibits  the  same  shortcomings. 

The  second  approach  is  to  search  in  some  manner  for  the 
coefficients  of  a  set  of  differential  equations  of  specified 
order,  the  response  of  which  is  sufficiently  close  to  that 
of  the  original  system  when  both  are  driven  by  the  same 
inputs.  Sinha  and  Pille  [10]  proposed  a  reduction  technique 
based  on  the  iterative  application  of  the  matrix  pseudo 
inverse  algorithm  [11]  to  determine  a  model  of  specified 
order  which  minimizes  the  sum  of  the  squares  of  the  errors 
between  the  responses  of  the  original  system  and  the  reduced 
order  model  to  a  given  input.  The  main  drawback  of  this 
method  is  that  the  objective  function  to  be  minimized  is 
restricted  to  be  the  sum  of  the  squares  of  the  errors.  Sinha 
and  Bereznai  [12]  presented  a  method  which  minimizes  a 
specified  error  criterion  for  a  given  reduced  order  model 
of  the  original  system,  based  on  the  pattern-search  algorithm 
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of  Hooke  and  Jeeves  [13].  Although  this  method  provides 
more  flexibility  than  that  of  Sinha  and  Pille,  it  generally 
requires  considerably  more  computational  time  due  to  the 
poor  convergence  properties  of  the  pattern-search  algorithm. 

The  third  category  involves  application  of  the  theory 
of  continued  fractions.  Methods  involving  this  approach 
have  been  developed  by  Chen  and  Shieh  [14], 

Sinha  and  DeBruin  [15]  and  Fellows  et  al  [16]  have 
established  the  fact  that  among  the  methods  previously 
mentioned,  the  approach  by  continued  fraction  expansion  is 
generally  the  best  for  linear  model  simplification. 

A.  SIMPLIFYING  A  TRANSFER  FUNCTION 

The  general  nature  of  a  control  system  is  that  of  a  low 
pass  filter.  Therefore,  model  simplification  should  con¬ 
centrate  on  the  steady  state  aspects  of  the  response  with 
the  transient  portion  given  secondary  consideration.  As 
previously  shown  in  Chapter  II,  the  Cauer  Second  Form 
exactly  characterizes  these  miens. 

Given  the  nth  order  original  system  transfer  function: 

n-1  t 

L  bis1 

T(s)  =  — — - - -  ,  (3-1) 

sn  +  I  a.s^ 

j=0  ^ 

where  an  mth  order  simplified  model  of  the  system  (where 
m  is  strictly  less  than  n)  is  desired,  the  polynomials  in 
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equation  (3-1)  are  rewritten  in  ascending  order; 


T(s)  = 


bQ^b^s^.  •  •  • 


K  ^n-2J.K  ^-1 

.  .  . D  „S  +D  ,  S 

x-2  n-i 


ao+ais+* 


x-2  ,  n-1,  n 

•'ax-2s  +an-ls  *3 


(3-2) 


and  expanded  into  a  continued  fraction: 


T(s)  = 


hl  + 


h2  + 


h3  + 


h2n-i  +  - § - ‘ 

h2n  (3-3) 

An  mth  order  simplified  model  is  obtained  by  keeping  the 
first  2m  quotients  of  the  expansion,  omitting  the  remainder; 

T(s)  =  _ 1 _ 

h-^  +  _ s _ 

hg  +  _ s _ 


h2m-l  + 


2m 


(3-4) 


and  performing  the  inversion  of  the  truncated  fraction.  The 
inversion  technique  presented  in  Chapter  II  can  be  used  for 
this  purpose. 
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Consider  the  seventh-order  system,  representing  the 
control  system  of  the  pitch  rate  of  a  supersonic  transport 
aircraft  [10  ],  described  by  its  transfer  function: 


T(s)  = 


375000(s+.  08333) 


s7+83.635s6+4097s5+70342s4+853703s3 


+  2814271s  +  3310875S+281250 


By  continued  fraction  expansion: 


T(s)  = 


9.00036  + 


-.486286  + 


(3-5) 


-.036856  + 


78496.2032+  _ s 

.00071478 


(3-6) 


Suppose  a  second-order  simplified  model  is  desired.  Equa¬ 
tion  (3-6)  has  fourteen  quotients.  For  the  desired  system, 
the  first  four  quotients  are  kept  with  all  others  discarded. 
The  truncated  continued  fraction  becomes : 

T1(s)  =  _ 1 _ 

9.00036  +  _ s _ 

-.  486286  +  _ s _  (3-7) 

-.036856  +  s 


.616185 
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and  converted  into  transfer  function  form; 


Tx(s) 


.1299S+. 01105 
s2+1.14644s+. 09941 


C  3-8 ) 


The  block  diagrams  of  equations  (3-5)  and  (3-8)  in  the 
Cauer  Second  Form  are  shown  in  Figures  3.1  and  3.2,  respec¬ 
tively.  The  unit  step  and  impulse  responses  of  the  original 
and  simplified  systems  an  compared  and  shown  in  Figures 
3 . 3  and  3.4. 


B.  STATE  EQUATION  SIMPLIFICATION 

The  method  of  system  simplification  just  presented  is 
especially  advantageous  when  converted  into  state  space 
form.  In  Figure  2.3,  a  name  for  each  state  variable  is 
given  after  each  integrator,  shown  in  Figure  3.5,  from 
which  the  state  equations  and  output  equation  can  be 


directly  written. 


h4hl  h6hl  "•h2nhl 

VW  VW  ••*h2n(hl+h3) 

h4(hl+h3)  h6(hl+h3+h5) ’ ’ *h2n(hl+h3+h5) 
•  •  • 

•  •  • 

•  •  • 

hl|(h1+h3 )  hg  (h1+h3+h5  )..  'h2n^hl+h3+ ' 
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486 


Figure  3.1.  Block  Diagram  of  7th  Order  System  from  Continued  Fraction 
Expansion 


98  <7 
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Figure  3.2.  Block  Diagram  of  2nd  Order  System  from  Truncated  Continued 
Fraction  Expansion 


Figure  3.3.  Unit  Impulse  Response  Comparison  of 
7th  and  2nd  Order  Systems 


Figure  3.4,  Unit  Step  Response  Comparison  of 
7th  and  2nd  Order  Systems 


State  Variables 


1 


1 


1 


D-  = 
~P 


1 


1 


and  =  h  z  *  , 
p  ~p  ~p’ 


c  =  [h„  h 
P  2 


h2m] 


(3-14) 


(3-15) 


(3-16) 


As  an  example,  consider  the  seventh  order  system  described 
by  the  transfer  function: 


C(s)  _  1441.53s3+78319s^+525286.125s+607693.25 
R<'S)  s7+112.04s6  +  375  5. 9  2s5  +  39  736 . 73s4 


(3-17) 


+  363650. 56s3  +  759894.19s2  +  683656.25s  +  617497. 375 


Arranging  the  polynomials  in  ascending  order  and  expanding 
into  the  continued  fraction  yields: 


0.2516 


Figure  3.7.  7th  Order  System  Reduction  to  2nd  Order 


Order  Systems 


C(s) 

rTsT  ' 


1 


1.016133  +  _ s 

4.054112  +  _ s _ 

—.067134  +  s 


(3-18) 


595660.646  +  s _ 

.0000757 


and  equations  (3-9)  and  (3-10)  are  formulated  from  the 
quotients.  A  simplified  model  of  second  order  is  desired, 
therefore,  the  state  and  output  equations  are  partitioned 
as  indicated  in  Figure  3.7.  The  simplified  transfer 
function  is: 


Cp(  s ) 
rpTsT 


. 250367S+1. 035264 
s2+.509768s+l. 051966 


( 3-19) 


Unit  step  and  impulse  responses  of  the  original  seventh  and 
simplified  second  order  systems  are  shown  in  Figures  3.8  and 
3.9  respectively.  It  is  observed  that  the  results  of  ex¬ 
pressing  the  seventh-order  system  by  a  second  order  model 
are  satisfactory. 

It  should  be  pointed  out  that  a  stable  transfer  func¬ 
tion  may  produce  an  unstable  simplified  function  because  the 
method  of  continued  fraction  expansion  approximation  does 
not  necessarily  guarantee  a  stable  system. 
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IV.  DESIGN  OF  OPTIMAL  LINEAR  CONTROL  SYSTEMS 
WITH  PRESCRIBED  EIGENVALUES 

A.  INTRODUCTION 

Consider  the  control  of  a  plant  with  dynamics  described 
by  a  set  of  first  order,  time-invariant  linear  differential 
equations  of  the  form 

x  =  Ax  +  Bu,  (4-1) 

where  x  is  the  nth-order  state  vector,  A  is  the  (nxn) 
plant  matrix,  u  is  the  scaler  control  and  B  is  the  (nxl) 
input  matrix.  The  output  is  defined  as 

y  =  Cx  ,  (4-2) 

where  C  is  the  (lxn)  output  matrix. 

A  linear  feedback  control  law  is  assumed,  and  of  the  form 
u  =  G*x . +  (4-3) 

There  are  mainly  two  separate  approaches  in  the  deter¬ 
mination  of  the  feedback  control  matrix,  G*,  corresponding 
to  the  system  under  consideration;  1  -  optimal  control  and 
2  -  modal. 

In  the  optimal  control  approach  a  performance  index  is 

considered  which  is  to  be  minimized  in  the  design  of  a 

All  states  are  available  or  an  observer  or  Kalman  filter 
is  used  to  obtain  the  unknown  states. 


system.  Assuming  a  performance  index  can  be  defined  that 
represents  most  of  the  design  requirements,  "the  solution 
to  the  optimal  control  problem  can  usually  be  obtained  only 
by  numerical  methods  that  yield  solutions  to  only  a  parti¬ 
cular  problem"  [17].  If  solutions  are  sought  to  more  than 
one  numerical  problem,  simple  performance  indices  must  be 
defined,  which  often  do  not  satisfy  many  of  the  design 
requirements.  Therefore,  the  choice  of  a  performance  index 
must  fall  somewhere  between  a  realistic  criterion  and  one 
that  is  mathematically  tractable. 

A  quadratic  performance  index  will  be  considered  as  a 
criterion  for  designing  linear  systems,  of  the  form 

CO 

J  =  h  f  [xT  Q  X  +  Ru2]dt  ,  (4-4 

o 

where  Q  is  a  diagonal  non-negative  definite  (nxn)  matrix, 
and  R  is  a  positive  scalar. 

In  the  modal  approach,  G*  is  chosen  so  that  the  closed 
loop  system  achieves  the  prescribed  eigenvalues.  Equations 
(1)  and  (3)  together  yield 

x  =  (A  +  BG* ) x  . 

If  Q  and  R  are  given  in  the  optimal  control  approach, 
then  the  eigenvalues  of  the  closed  loop  system  are  uniquely 
determined,  which  may  not  realize  the  required  performance 
characteristics  or  desired  degree  of  stability  for  the 
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system.  Using  the  modal  control  approach,  a  feedback  control 
matrix  can  be  found  that  will  give  the  system  the  desired 
eigenvalues.  This  matrix  is  usually  not  unique,  and  it  is 
not  possible  in  a  practical  sense  to  find  one  that  is 
"better"  than  it's  predecessors,  since  a  performance  measure 
is  generally  not  known  that  corresponds  to  a  given  feedback 
control  matrix.  Therefore,  it  is  necessary  to  find  a  method 
for  determining  the  matrix,  G* ,  that  simultaneously  satisfies 
the  desired  eigenvalues  and  minimizes  a  given  performance 
index. 

In  addressing  this  problem,  Chang  [18]  and  Tyler  and 
Tuteric  [19]  have  applied  the  root  locus  method  to  single¬ 
input,  single-output  and  multivariable  systems,  respectively. 
The  method  lacks  a  rational  computational  procedure  for 
determining  the  elements  of  the  weighting  matrix,  Q,  to  meet 
a  set  of  prescribed  closed  loop  eigenvalues.  Anderson  and 
Moore  [20]  presented  a  restrictive  method  whereby  a  set  of 
eigenvalues  may  be  located  to  the  left  of  a  line  parallel 
to  the  imaginary  axis  in  the  complex  plane.  Chen  and  Shieh 
[14]  presented  a  method  using  sensitivity  analysis. 

Solheim's  [17]  method  of  a  diagonalized  (decoupled)  system 
becomes  complicated  when  the  system  contains  either  complex 
or  multiple  eigenvalues.  Systems  that  cannot  be  diagonalized 
add  further  to  the  complication  of  the  method. 

The  method  developed  here  takes  advantage  of  the 
properties  of  the  Cauer  Second  Form,  is  approached  in  a 
simplistic  manner,  and  is  easy  to  implement  computationally. 


B.  TRANSFORMATION  TO  PHASE  VARIABLE  FORM 


Consider  an  nth  order  linear  system  of  the  form 
e  =  Se  +  Tf  ,  (4-5) 

with  output 

d  =  We  .  (4-6) 

Silverman  [21],  et  al.,  have  shown  that  if  the  system  is 
controllable,  then  there  exists  a  non- singular  transformation 
matrix  which  takes  an  arbitrary  state  variable  system  to 
phase  variable  (canonical)  form.  (see  Appendix  A) 

x  =  Ax  +  Bu  (4-1) 

y  =  C  x  ,  (4-2) 


where 


If  the  system  is  not  controllable,  the  phase  variable 
(canonical)  form  may  still  be  obtained  from  the  system 
transfer  function 


Y(s)  _ 
uTsT 


1  8  Si_1 

i  =  l  BiS _ 

sn  +  Z  cl  .  s1"1 
i  =  l  nl 


(4-7) 


Once  the  system  is  in  phase  variable  form,  Chen  and  Shieh 
[22]  have  shown  that  the  equivalent  system  in  Cauer  Second 
Form  is  easily  written  as 

z  =  Hz  +  Vu  (4-8) 

y  =  C*z  ,  (4-9) 

where  the  two  forms  are  related  by  a  linear  nonsingular  trans¬ 
formation  matrix,  P, 

z  =  Px  .  (4-10) 

The  matrix  P  is  an  (nxn)  upper  triangular  matrix. 
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The  performance  measure  under  consideration  becomes 

QO 

J  =  \  /  Cz^Qz  +  u^RuJdt  (4-11) 

o 

where 

Q  =  (P_1)T  Q  P"1  . 

From  optimal  control  theory,  the  Hamiltonian 

H  =  isC2TQz  +  uTRu]  +  PT[Hz  +  Vu]  ,  (4-12) 


where  P  is  the  set  of  Lagrange  Multipliers  (also  called  the 
costate  or  adjoint  vector) .  For  the  Hamiltonian  to  be 
globally  minimized,  assuming  no  bounds  on  admissible  states 
and  control  values,  it  is  necessary  that  3H/3u  =  0  and 
32H/9u2>0. 

3  H/3  u  =  Ru  +  VTp  =  0  (4-13) 

implies 

u*  =  -R  lVTp  ,  and  (4-14) 

32H/3u2  =  R>0  ,  (4-15) 


since  R  was  defined  as  a  positive  scalar.  Included  in  the 
necessary  conditions  for  optimality  are 


z  =  Hz  +  Vu*  (4-16) 

X  -  •p/. 

3 H/3  z  =  -P  =  Qz  +  H  p  .  (4-17) 
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Combining  equations  (4-14),  (4-16)  and  (4-17)  yields  a  set  of 
2n  linear  differential  equations  forming  the  canonical  sys¬ 
tem  in  Cauer  Second  Form. 


1 

»  N  • 

l _ 

vr"1vt 

1 - 

Nl 

1 _ 

m 

A 

P 

II 

-Q 
_  ~ 

-ht 

- 1 

<0.1 

_ 1 

(4-18) 


It  remains  to  be  shown  that  this  form  is  useful  in  obtaining 
optimal  closed  loop  solutions  that  correspond  to  a  set  of 
prescribed  eigenvalues. 


C.  SIMILAR  EIGENVALUES 


Consider  the  2nth  order  cononical  system  in  phase 
variable  form: 


A 

■N. 

-Q 


BR'1 


(4-19) 


It  is  known  that  this  system  possesses  n  eigenvalues  with 
negative  real  parts  and  n  eignevalues  with  positive  real 
parte,  and  that  they  are  located  symmetrically  about  the 
imaginary  axis  in  the  complex  plane  [  17  ].  The  eigenvalues 
of  the  optimal  feedback  system 


x  =  (A  +  BG*)x 


(4-20) 


are  identical  to  those  with  negative  real  parts  in  the 
canonical  system.  It  is  possible,  therefore,  to  focus 
on  the  canonical  system  where  the  dependence  eigenvalues 
on  the  weighting  matrices  Q  and  R  can  be  studied  without 
solving  the  matrix  Riccati  equation.  The  problem  is  in 
determining  a  weighting  matrix,  Q,  such  that  the  system 
attains  the  prescribed  set  of  eigenvalues. 

The  canonical  system  in  Cauer  Second  Form  can  be  ob¬ 
tained  from  the  phase  variable  form  using  the  linear  trans¬ 
formation 

z  =  P  x  ,  (4-10) 


where  P  is  a  nonsingular  matrix.  We  have 


A  -BR'XB  ~ 

1 

Xi 

1 _ 

-Q  -AT 

- 1 

0.1 

J 

- 1 

I  N  • 

_ 1 

P  0 

.  1  fi 

A  -BR  B 
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_  P_ 

II 

_0  (P_1)T_ 

T 

L-Q  -A1  J 
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(4-19) 


(4-20) 


(4-21) 


where  each  sub-matrix  of  M,  N,  and  P  are  known  to  be  Cnxn) 
square  matrices.  It  is  easily  seen  that  P  is  non-singular, 


and  that 


(P_1)T 


P  P_1  0 


(P-W 


I 


(4-25) 


Therefore , 


P  M  P"1  =  N 


(4-26) 


shows  the  similarity  of  the  M  and  N  matrices.  Two  similar 
matrices  have  the  following  properties: 


1.  Their  determinants  are  the  same. 


det  M  =  det  N 


(4-27) 


2.  Their  traces  are  the  same. 


Tr  CM]  =  TrCN]  (4-28) 

k  k 

I  m . .  =  I  n..  (4-29) 

i  =  l  11  j=l  133 

3.  Their  characteristic  equations  are  the  same. 

detCXI-M]  =  detCXI-N]  =  0  ,  (4-30) 

where  X  is  an  arbitrary  variable.  Since  their  characteristic 
equations  are  the  same,  the  eigenvalues  of  M  and  N  must  be 
identical.  It  is  now  known  that  the  Cauer  Second  Form  and 
phase  variable  system  matrices  are  similar  in  that  they 
possess  identical  eigenvalues. 

D.  DEVELOPMENT  OF  THE  PRESCRIBED  EIGENVALUE  PROBLEM 

Initially  given  is  a  known  linear  system  with  dynamics 
described  by  either  a  set  of  first  order  differential  equa¬ 
tions  or  its  transfer  function.  It  is  desired  to  find  the 
optimal  feedback  control,  u* ,  such  that  the  performance 
measure 


uT[l]u]  }  dt 


(4-31) 


is  minimized,  where  the  eigenvalues  of  the  optimal  system 
are  specified  as 


1.  Evaluation  of  the  State  (H)  and  Linear 
Transformation  (P)  Matrices 

Assume  the  transfer  function  of  the  known  system  is 


given 


T(s) 


n 

E  B  ,s 
i=l  1 


i-1 


„  n 
sn  +  E 
i=l 


i-1 


(4-32) 


The  system  matrices 


in  phase  variable  form  are: 


(4-33) 


C  —  C8  -i  8^  •  •  •  •  8  3 
-  l  i  n 


(4-33) 


From  equation  (4-32),  the  Routh  array  is  formed: 
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d2  . 
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(4-34) 


In  matrix  notation,  the  Routh  array  becomes  [r^] 


where 


ieCl ,  2 ,  2 ( n+1) ] 

j  e[l,  2,  n+1]  , 

and  elements 

r2( n-k) +3 ,  k  =  1 

r2(n-k)-4 ,  k  =  0 

kc[l ,  2 ,  .  .  .  .  ,  n+l]  . 


In  general,  the  (2L+l)th  row  of  the  Routh  array  is  the  Lth 
row  of  P. 
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p  = 


31 


32 

r33  - 

3  ,n 

51 

52  . 

5  ,n-l 

. 7 ,n-2 

2n+l , 2 
2n+l ,  1 


(4-35) 


P  being  an  (nxn)  upper  triangular  matrix,  will  always  have 
an  inverse,  P-1,  which  can  be  quickly  and  efficiently  deter¬ 
mined  by  digital  computer.  The  H  matrix  formulation  becomes: 


H  =  P  A  P 


-1 


(4-36) 
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h33  . 
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(4-37) 


The  elements  of  the  H  matrix  can  also  be  obtained  more 
easily  and  directly  from  the  first  column  of  the  Routh  array. 
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Let 


h. 

i 


U 


i+1,1 


(4-38) 


The  h^'s  correspond  to  the  quotients  of  the  continued  frac¬ 
tion  expansion  of  the  transfer  function  in  (32). 


T(s)  = 


hl  + 


h2  + 


h3  + 


h4  +  . 


(4-39) 


The  H  matrix  then  becomes  as  shown  in  Figure  4.1.  The 
regular  pattern  of  the  elements  enable  the  H  matrix  to  be 
obtained  by  inspection  once  the  fu  ' s  have  been  determined 
from  either  (4-34)  and  (4-38),  or  (4-39). 

The  matrices  V,  and  C*,  are  easily  obtained: 


V  =  P  B  = 


1 


(4-40) 


Figure  4.1.  State  Matrix,  H,  in  Cauer  Second  Form 


(4-41) 


r 


c*  =  C  p'1  =  [h2  h4  hg  -  h2n] 

B.  Evaluation  of  Q 


Q  =  (P_1)T  Q  P"1 
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(4-42) 


The  canonical  system  in  Cauer  Second  Form  (4-23)  has  now 
been  obtained,  with  the  numerical  values  of  the  elements 
of  Q  still  to  be  found.  This  system  will  be  compared  to 
a  non-optimized  system  with  the  prescribed  eigenvalues  also 
in  Cauer  Second  Form.  The  desired  system  in  phase  variable 
form  is : 


x  =  A*x  +  B*u 
y  =  E  x 

Formulation  into  Cauer  Second  Form  yields 
x  =  H*z  +  V*u 
y  =  E*z 


(4-43) 

(4-44) 


(4-45) 

(4-46) 
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In  a  "nonoptimized"  system  (Q  and  R  set  equal  to  0),  the 
canonical  system  appears  as: 


z 

“H* 

0 

« 

_  p_ 

_0 

-(H*)T  _ 

where 


(4-47) 


“H* 

0 


(4-48) 


possesses  the  n  prescribed  eigenvalues  with  negative  real 
parts  and  n  eigenvalues  with  positive  real  parts,  symmetric 
about  the  imaginary  axis  in  the  complex  plane . 

By  equating  the  characteristic  polynomials  of  N  and  N* , 
(4-23)  and  (4-47)  respectively,  it  is  now  possible  to 
determine  the  elements  of  the  Q  weighting  matrix. 


det[sI-N]  =  det[sI-Nft] 


det 


sI-H 


Q 


vr"1vt 

sI+HT 


(4-49) 


det 


sI-H* 

0 


0 

sI+(H*)^ 


(4-SO) 
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E.  DETERMINATION  OF  THE  OPTIMAL  CONTROL  LAW 


Finding  the  elements  of  the  weighting  matrix,  Q,  is 
obviously  a  non-trivial  matter  for  all  but  the  lowest  order 
systems.  The  method  developed  for  obtaining  these  values 
is  based  upon  a  succession  of  matrix  building  blocks,  which 
are  individually  computationally  simple. 

Starting  with  the  matrix,  H, 


IK 

II 

ch..] 

(4-51) 

define 

a  new  matrix,  T, 

T  = 

(4-52) 

where 

tij 

=  hi3  -  h.. 

i 

*  j 

(4-53) 

and 

n 

tii 

=  E  t .  . 

j-i+1  13 

i 

i  n 

(4-54) 

The  matrix,  T,  therefore,  is  an  (nxn)  upper  right  triangular 
matrix,  where  the  diagonal  elements  are  equal  to  the  sum 
of  all  other  elements  in  the  same  row.  The  next  "building 
block"  matrix,  G,  is  defined  by: 

G  =  [gi-]  (4-55) 

gijX  =  1.0  (4-56) 
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for  ietl,  2 ,  . . . ,  n]  , 


’i»2  = 


for  jctl,  2,  n-1]  , 


n-j  +  2 

S-i-i  =  E  v  x  Sv  -i-i5 

k=i+l  1,K  K,:!  1 


(4-57) 


(4-58) 


for  je[3,  4,...,  x]  and  iE[l,  2,...,  n-j+1],  where  the 
index  j  is  held  fixed  for  each  summation  over  the  index  i. 
The  matrix,  G,  is  an  (nxn)  upper  left  triangular  matrix 
characterized  by  the  first  column  being  all  ones.  One  more 
matrix  needs  to  be  defined  at  this  point.  Let  the  matrix, 
W,  be  defined  as: 


(4-59) 


where 


*  (-1>:  x  *i,j 


(4-60) 


for  i,  j  and  ke[l,  2,  ...,  n] 


The  matrix,  W,  is  therefore  a  tridimensional  array  with  each 
"level"  an  upper  left  triangular  matrix.  Examples  of  the 
matrices,  T,  G  and  W  are  shown  below.  Although  not  evident 
at  this  point,  the  G  and  W  matrices  will  be  used  heavily  in 
obtaining  the  values  of  the  elements  of  the  Q  matrix. 


From  the  linear  transformation 


Q  =  (P~1)T  Q  PT  ,  (4-61) 

it  is  observed  that  once  the  element  is  known,  the 

remaining  elements  in  the  same  row  (column)  can  then  be 
obtained  through  a  process  of 


(4-62) 
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linear  combinations  of  previously  determined  values.  The 
Q  matrix  is  thus  found  in  the  following  manner: 


+  for  n=l ,  Qn  =  (H2H1)2  -  <h2hx)2. 
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(4-67) 
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for  i^j  or  j/n  , 


~  n-1 

Q .  =  -  E  Q .  . 

in  j=l  ^ 


(4-67) 


and 
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E  E  h. .(h. .-h, .)]  (4-69) 

L=1  j=l  1]  33  13 


The  values  of  the  diagonal  elements  of  Q  (other  than  Q^) 
generally  involve  varying  linear  combinations  of  already 
determined  values,  more  easily  expressed  in  terms  of  the 
and  W  matrices,  rather  than  the  P  and  H  matrices.  To  aid 
in  determining  the  values  along  the  diagonal  the  following 
labels  are  provided  for  G  and  W. 
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=  P 
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CD  I 


These  two  arrays  will  provide  the  coefficients  of  each  "Q^x 
s  "  necessary.  Assume,  for  example,  that  a  5th  order 
system  is  being  considered  in  equation  (4-50).  The  results 
of  the  expansion  of  both  sides  of  the  equation  results  in: 


«  2i  11  2i 

a  +  E  a.s  =  b  +  E  b-s  (4-72) 

°  i=l  1  °  L=1  1 

where  an  =  1,  b  =  1.  By  equating  a^_^  and  b^_1  for 
ie[l,2,  ...»  n],  it  is  possible  to  obtain  Q... 


i  .e . 


(4-73) 


a  and  b  are  the  coefficients  of  s°.  r'rom  (4-70),  it  is 
o  o 

observed  that  the  only  element  in  the  s°  column  that  is 
non-zero  is  g15»  which  appears  in  row  1.  The  1  "indicates" 
that  it  is  necessary  to  only  look  in  "level"  1  of  W,  in  the 
column  corresponding  to  s°.  This  yields  only  the  single 
element  Therefore,  the  coefficient  of  will  be 

Sis  x  W1S1  •  (4-7 


What  remains  to  be  determined  are  the  coefficients  of  s71 
not  involving  the  Q^’s.  Because  of  the  symmetry  of  the 
eigenvalues  of  both  the  system  being  "optimized"  and  the 
system  with  the  prescribed  eigenvalues,  the  remaining 
coefficients  (those  not  involving  a  are  easily 

determined  from: 
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det[s!-H]  x  detCsI  HT] 


(4-75) 


and 


det[sI-H*]  x  det[sI+(H*)T] , 


(4-76) 


which  give: 


(  I  a-  s1)(  Z  (-l)ka.  s1)  (4-77) 

i=0  i=0  1 

and 


(  Z  a • *s1) (  Z  (-l)ka.*s1)  (4-78) 

i=0  1  i=0  1 


where  k-i+1  for  n  even  and  k  =  i  for  n  odd.  The  indicated 
multiplication  in  (4-77)  and  (4-78)  result  in: 

Z  (-l)ka-s21  (4-79) 

i=0  1 

and 

Z  (-l)ka.*s21  (4-80) 

i  =  0  1 


respectively,  where  the  same  conditions  are  imposed  on  "k" . 

Returning  to  the  5th  order  example,  the  coefficients  of 
s°  are  -aQ  and  -aQ*.  To  obtain  <5)^  is  now  a  matter  of 
solving  the  equation: 
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(4-82) 


On  ~  Or 


'll 


Ti 


15 


g151) 


With  known,  it  is  a  simple  matter  to  obtain  the 
remaining  elements  in  the  first  row  (and  column)  using 
(4-67)  and  (4-68) . 

~  2  2 
To  find  Q22’  all  the  coefficients  of  s  .  s  results 

in  three  separate  ways : 


(1) 

(2) 

(3) 


0 

s 

1 


X  s 

X  s 

X  s 


2 

1 

0 


(4-83) 


The  multiplicand  indicates  which  columns  of  G,  (4-70),  is 
of  interest.  Any  non-zero  element  in  G  tells  which  level 
of  W  is  of  interest.  The  multiplier  is  the  indicator  for 

the  column  of  interest  in  the  array,  W.  Therefore,  for 

0  2 
s  x  s  : 

g15  X  (W131  Q11  +  W231  Q12  +  W3  31  Q135  (4 

for  s1  x  s^: 


gl,4  X 

(W141 

5u 

+  W241 

Ql2)  - 

g2 , 4  X 

(W14  2 

*21 

+  W242 

Q22) 

(4-85) 

80 


--■n* 


_  2  0 
for  s  x  s  : 


gl , 3  X  (W151  Qll)  + 


g2 , 3  X  (W152  Q21)  + 


g3 , 3  X  (W153  Q31} 


(4-86) 


Obtaining  Q22  is  now  a  matter  of  solving  the  equation  given 
by  the  combination  of  (4-79),  (4-80),  (-84),  (4-85)  and 
(4-86) : 


3  2  2 

gl5  ^  ^  Qi-?)  +  ^  2  ii  £  W-ju-;  Q^,*) 


i  =  l 


i31  yli' 


l=i  i=1  "i4j  vji' 


(}Zsl  gj,3  W15j  Qji)  *  °1  =  "°1*  » 


(4-87) 


for  Q22. 

Once  Q22  has  been  found,  the  remainder  of  the  2nd  row 
(and  column)  can  be  obtained  using  (4-67)  and  (4-68). 

In  a  like  manner,  the  respective  equations  to  be  solved 
for  Q33,  Q,4  and  Q55  are: 


gi>5  V  +  i=1Wi21  V 

*  (3i  gi.3  j/uj  V  *  2  j^iaj  5ji> 


5 

Z 

j=l 


*  (-E  gj.l  J^iSj  -  °2  *  -°2* 


(4-88) 
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*  Sj,3  J/ilj  V  +  <.^,2  V 


*  < Jj.  *3,1  j^isj  V  *  °3  =  “a8  • 


(4-89) 


(j=lgj)1  i=lWilj  Qji>  "  °4  =  ”°4*  ' 


(4-90) 


The  underscore  in  (4-88),  (4-89)  and  (4-90)  indicates  the 
term  that  contains  Q33)  0442  and  Q55  »  respectively . 

The  entire  Q  matrix  is  now  known.  By  the  inverse 
transformation  of  (4-42), 

Q  =  PT  Q  P  .  (4-91) 

Q  along  with  A,  B  and  C  of  (4-1)  and  (4-2)  are  used  to 
obtain  the  matrix  Riccati  equation  steady  state  solution. 

0  =  KA  +  ATK  -  KBR-1BTK  +  Q  .  (4-92) 

Once  K  has  been  determined,  the  optimal  control  law  is 
given  by 

u*  =  -R~1BKx  .+  (4-93) 

Once  u*  has  been  determined,  an  inverse  non-singular  trans¬ 
formation  can  be  performed  to  take  the  phase  variable  form 

+  It  is  now  known  that  G=-R~1BTK  in  (4-20). 


» 


i 
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back  into  state  variable  form.  (See  Appendix  A) 


F.  ILLUSTRATE  EXAMPLES 

1 .  Odd  Order  System  (Third) 
Given  the  system 


X1 

0 

1 

— 

0 

x2 

- 

0 

0 

i 

_*3_ 

-13 

-19 

-7 

y  =  Cl  0 

0] 

rxii 

> 

X2 


x. 

0 

1 

X2 

+ 

0 

1 

X 

Cl) 

1 _ 

1 

(4-95) 


find  the  optimal  control  law,  u*,  such  that  the  quadratic 
performance  index  (4-31)  is  minimized,  where  the  eigenvalues 
of  the  system  are  specified  as 


s 


1 


-3, 


=  -*+, 


-6. 
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Forming  the  Routh  array  yields 

13  19  7  1 

10  0  0 
19  7  1 

-7/19  -1/19  0 

30/7  1 

1/30  0 

1 
0 

from  which  the  third,  fifth  and  seventh  row  are  extracted  to 
form  P. 

“19  7  1 

P  =  0  30/7  1 

_0  0  1 

and 


(4-98) 


From  (4-97) 

H=PAP~1= 


and  (4-98), 
-13/19 
-13/19 
-13/19 


H  and  V  are 
637/570 
-63/19 
-63/19 


calculated 


(4-99) 


84 


1 


(4-100) 


The  set  of  prescribed  eigenvalues  yields  the  system: 


(4-101) 


From  the  Routh  array 


72 

54 

13 

1 

0 

0 

54 

13 

1 

-13/54 

-1/54 

0 

115/13 

1 

1/115 

0 

1 

0 


(4-102) 


(4-103) 
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extracting  rows  3,  5  and  7  yields  the  P*  matrix: 


~  54 

13 

1“ 

"  1/54 

-169/6210 

1/115  ~ 

II 

* 

a*  > 

0 

115/13 

1 

,  (P*)"1  = 

0 

13/115 

-13/115 

_  0 

0 

1  _ 

_  0 

0 

1 

(4-104) 


from  which  H*  and  V*  are  calculated, 

-»*  -X 


“  -4/3 

676/345 

-72/115 

H*=p*A*(P*)~1  = 

-4/3 

-2574/621 

1980/1495 

-4/3 

-2574/621 

-11254/1495 

(4-105) 


V*  s  p*B* 


1 


1 


1 


(4-106) 


Formulation  of  equation  (4-50)  yields 


si  -  H 

vr_1vt 

sI-H* 

0 

det 

Q 

sI+HT 

=  det 

0 

T 

sl+ff 

86 


det 


=det 


s+13/19 

-637/570 

13/30 

1 

1 

1 

13/19 

s+63/19 

-9/7 

1 

1 

1 

13/19 

63/19 

s  +  3 

1 

1 

1 

?11 

\l2 

\l3 

S-13/19 

-13/19 

-13/19 

^12 

q22 

Q23 

637/570 

s-63/19 

-63/19 

?13 

Q23 

Q33 

-13/30 

9/7 

s-3 

- 

"s  +  4/3 

-676/345 

72/115 

0 

0 

0 

4/3 

s+2574/621 

-1980/1495 

0 

0 

0 

4/3 

2574/621 

s+11245/1495 

0 

0 

0 

0 

0 

0 

s-4/3 

-4/3 

-4/3 

0 

0 

0 

676/345 

9-2574/ 

-2574/ 

621 

621 

0 

0 

0 

-72/115 

1980/  s 

-11245, 

1495 


1495 
(4-107) 


It  is  desired  to  determine  the  values  of  Q.  Brute  force 
enumeration  of  the  determinants  and  equating  coefficients 
of  like  powers  of  s  would  eventually  lead  to  the  solution. 
Using  instead  the  method  developed,  from  H(not  H* )  evolve 
the  T,  G,  and  W  matrices. 


T  = 


7.0000  4.4333  2.5667 
0  4.2857  4.2857 
0  0  0 


(4-108) 
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G  = 


1  7.0000 

1  4.2857 

1  0 


19.0000 

0 

0 


(4-109) 


and  W  = 


-1  7.0000 
-1  4.2857 
-1  0 


-19.0000 

0 

0 


(4-110) 


Now  add  the  appropriate  labels  as  in  (4-70)  and  (4-71) 


for  G: 


1 

2 

3 


s 

1 

1 

1 


s  s 

7.0  19.0 

4.2857  0  ' 

0  0 


(4-111) 


for  W: 


9n 


912 

913 


■  1 
•1 
■1 


7.0  -19.0 

4.2857  0 

0  0 


(4-112) 


Still  necessary  are  the  coefficients  in  (4-79)  and  (4-80) 


n 


n 


(  I  a.s1)(  Z  (-1)  a.s1)  =  Z  (~l)Ko.s 
i=0  1  i=0  1  i=0  1 


2i 


(4-113) 


for  n  odd  (n=3),  k=i+l, 


(  Z  a.si)(  Z  (-1) 1+1a. s1)  =  Z  (-l)1+1a. s21  (4-114) 

i=0  1  i=0  1  1=0  1 


(s3+7s2+19s+13)(s3-7s2+19s-13) 


=  s6  -  11s4  +  179s2  -  169, 


(4-115) 


and 


3.3. 

(  Z  ct  •  *s1)  (  Z  (-l)1+1a.*s1) 
i=0  1  i=0  1 


3 

Z  (-1) 
i  =  0 


l+l 


a . 


2i 


(4-116) 


(s3+13s2+54s+72)(s3-13s2+54s-72) 
=  s6  -  61s4  +  1044s2  -  5184  . 


(4-117) 


Now,  by  equating  coefficients  of  s3,  it  is  possible  to 
obtain  Q^.  From  G,  the  only  non-zero  element  in  the  "s'"1 
column"  is  19,  which  corresponds  to  row  1,  and  therefore 
level  1  of  W.  From  level  1  of  W,  the  only  non-zero  element 
in  the  "s°  column"  is  -19,  with  row  label  Q^.  The 
coefficient  of  as  obtained  from  (4-81)  is  19(-19)  = 

-361.  The  solution  for  is  obtained  from  equations  (4-81) 

and  ( 4-82)  : 


'4 


1 
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-361  Qn  -  169  =  -5184 
Qn  =  (-5184  +  169)  /  -361 

=  13.892  .  (4-118) 


From  (4-65),  (4-67)  and  (4-68),  and  are  successively 

obtained: 


12 


12 

'22 


'll 


1 

"3077 


(13.892)  =  -22.69 


(4-119) 


Q13  S  "Qll  '  Q12  *  “13*892  “  (“22.69)  =  8.798.  (4-120) 


The  next  power  of  s  which  results  from  expansion  of  (4-107) 

2  2 
is  s  .  Equating  coefficients  of  s  ,  it  is  now  possible  to 

2  0  2  1  1 
obtain  Q22’  s  results  from  the  products  s  xs  ,  s  xs  , 

2  0 

and  s  xs  .  Therefore,  from  the  development  starting  at 
(4-83),  the  equation  to  be  solved  for  Q22  is: 


(j*x«i.3  J/ui  V  *  ‘ji/j.*  ^"121  V 


*  Su®*-1  V  *  °l  = 


(4-121) 


1 

1  \ 


i 


90 


Substituting  known  values,  equation  (4-121)  becomes 


S13(W111Q11  +  W2H  Q12  +  W311  + 


gl2(W121  <*11  +  W221  ^12)  +  g22(W122  ^21  +  W222  Q22) 


+  gll(W131  ^ll5  +  §21(W132  52i>  +  §31(W133  Q3l> 


+  °1  =  °1* 


(4-122) 


19. 0[ -1(13. 892)  -  K-22.69)  -  1(  8.793)] 


+  7.0  [7.0(13.892)  +  4 . 2 8 57 ( -2 2 . 69  )  ] 


+  4.2857  [ 7 . 0 ( -22 . 69 )  +  4.2857  Q22)]  +  1[ -19 ( 13 . 89 2  ) ] 


+  1  C  — 19(  —  22 .69)]  +  1[ -19 (8.798  )]  +  179  = 


1044 


Q22  =  84.155 


(4-123) 


From  (4-68) 


Q23  =  'Q21  '  ^22  =  _Q12  “  Q22 


-  (-22.69)  -  (84.155)  =  -61.465 


(4-124) 


a 
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4 

s  is  the  next  power  of  s  obtained  in  the  expansion  of  (4-107). 

4 

Equating  coefficients  of  s  ,  it  is  now  possible  to  obtain  the 

~  4 

equation  from  which  Q ^ 3  can  be  found,  s  results  from  multi- 
2  2 

plying  s  xs  .  Therefore,  from  G  and  W: 

*•» 

gnC-KQn)  -  1(Qi2>  -  HQ13>] 


+  g21C-l(Q21)  -  l(Q22)  -  l(Q23)3 

+  Ssi^^Qsi5  -  1(^32)  -  1(^33) 


-  a2  =  -a 2*  ,  (4-125) 

where  g^,  g21  and  g3i  equal  one.  Therefore,  (4-12  5) 

becomes 

3  3 

-  £  I  Q..  -  O  =  -o  *  (4-126) 

i  =  l  j  =  l  13  1  1 


Solving  (4-126)  for  Q33: 


<33 


•a2*  +  a2 


3  3  - 

Z  IQ.. 

i«l  j=l  13 

i+ j  /  6 


Since 

Q 


in 


n-1  - 

E  Qik 
k=l 


(4-127) 


(4-68) 
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0 


(4-128) 


n  ^ 

s  Qik  :  o  , 
k=l  1K 


and  (4-127)  simplifies  to 


«3k  =  *°2*  -  °2  -  «3k 


=  61-11-8.798*61.465  =  102.667 


Q  is  now  entirely  known. 


Q 


13.892 

-22.69 

8.798 


-22.69  8.798 

84.155  -61.465 

-61.465  102.667 


From  (4-91): 


Q 


5015 

0 

0 


0  0 

865  0 

0  50 


(4-129) 

(4-130) 


(4-131) 


(4-132) 


To  find  the  optimal  feedback  control  law,  u*,  it  is  necessary 
to  solve  the  matrix  Riccati  equation  (4-92),  where  Q,  A,  B 
and  C  are  as  given  in  (4-131),  (4-94)  and  (4-95).  The 
solution  yields: 
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u*  =  -C  59  35  6] 


or  u*  =  -59  -  35  Xj  -  6 

With  G*  as  given  in  (4-133),  x  =  (A  +  BG*)x 
becomes 


(4-133) 


(4-134) 


Comparing  (4-134)  with  (4-101),  it  is  observed  that  (A+BG*)= 
A*;  therefore,  the  desired  system  has  been  realized  with  the 
set  of  prescribed  eigenvalues. 

2 .  Even  Order  Example  (Fourth) 

Consider  the  linear  system  represented  by  the  transfer 
function : 


T4(s) 


_ s^9s+34 _ 

s4+12s3+48s2+80s+48 


(4.135) 


+ 


the  system  eigenvalues  are  -2, -2, -2, 


and  -6. 
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This  system  could  be  expressed  in  phase  variable  form  (4-31), 
thereby  obtaining  the  transformation  matrix,  P;  the  system 
can  be  realized  in  Cauer  Second  Form  (4-36).  Instead,  in 
this  example,  it  is  desired  to  obtain  the  quotients  which 
result  from  partial  fraction  expansion  of  the  transfer  func¬ 
tion.  In  Cauer  Second  Form,  T(s)  becomes: 


T(s)  = 


34  +  9s  +  s 

49+80s+48s2+12s3+s4 


hl  + 


V 


V 


h4  +  _ i 


h5  + 


h6  + 


h?  +  _ s_ 


(4-136) 
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where  the  quotients,  h^,  are: 


1.41176 

cn 

II 

23.07 

.505245 

h6  = 

.146914 

-4.6287 

tl 

r^ 

rC 

-281.808 

-.627918 

h8  = 

-.024241 

Substituting  these  values  into  the  matrix  in  Figure  4.1 
yields 


.71328 

-.88647 

.20741 

-.03422 

.71328 

2.01997 

-.47261 

.07798 

.71328 

2.01997 

2.91670 

-.48125 

.71328 

2.01997 

2.91670 

6.35004 

Again,  V  is  (and  in  all  cases  considered  will  be)  a  vector 
of  all  ones , 

(4-138) 

It  is  desired,  as  before,  to  find  the  optimal  control  law, 
u*,  such  that  the  quadratic  performance  index  (4-31)  is 
minimized,  and  that  the  optimal  system  realizes  a  set  of 
prescribed  eigenvalues.  The  eigenvalues  are  given  as: 
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-5 


S1 


-2, 


s 


2 


s3>  s4  =  “6  1  33 


(4 


With  the  zeroes  of  the  transfer  function  (4-136)  the  same, 
the  transfer  function  of  the  desired  system  with  the  pre¬ 
scribed  eigenvalues  becomes: 


T4*(s) 


s2+9s+34 

s4+19s3+138s2+435s+450 


(4 


By  continued  fraction  expansion  of  T*(s)  the  quotients 
obtained  are: 


H1  =  13.2353 
H2  =  .107635 
H3  =  -71.3206 
H4  =  -.088175 

from  which  H*  is: 


Hc  =  -1077.43 

D 

Hg  =  .004834 
H?  =  396.926 
Hg  =  -.024294  , 


H* 


1.42458 

1.42458 

1.42458 

1.4258 


-1.16702 

5.12168 

5.12168 

5.12168 


.06399 

-.28082 

-5.48975 

-5.48975 


-.32154 
1.41114 
27.58655 
17 . 94349 


(4 


-139) 


-140) 


-141) 
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Formulation  of  equation  (4-50)  yields  in  matrix  notation: 


sI-H 

l 

+ 

“ sI-H* 

0 

det  = 

Q 

T 

sl+n 

det 

0 

sI+(H*)T 

(4-142) 

which  is  to  be  solved  for  Q. 

As  before,  we  formulate  the  T,  G' and  W  matrices 
according  to  equations  (4-53)  and  (4-54);  (4-56),  (4-57) 
and  (4-58);  and  (4-60),  respectively. 


This  yields: 


o 

o 

C*J 

rH 

2.9064 

2.7093 

6.3500  “ 

0 

9.6614 

3 .3893 

6.2721 

0 

0 

6.8313 

6.8313 

_  0 

0 

0 

0 

”  1 

12.0000 

46 . 5882 

67.2941 

1 

9.6614 

23.1534 

0 

1 

6.8313 

0 

0 

1 

0 

0 

0 

(4-143) 


(4-144) 


the  notation  1  indicates  an  (nxn)  matrix  with  all  elements 
unity . 
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w  = 


sT 


-1 

12.0000 

-46.5882 

67.2941 

-1 

9.6614 

-23.1534 

0 

-1 

6.8313 

0 

0 

-1 

0 

0 

0 

(4-145) 


With  the 

proper  labels. 

(4-70)  and 

(4-71), 

for  G: 

3 

s 

2 

s 

s1 

0 

s 

1 

“  l 

12.0000 

46.5882 

67 . 2941“ 

2 

1 

9.6613 

23.1534 

0 

3 

1 

6.8313 

0 

0 

4 

_  1 

0 

0 

0 

(4-146) 


for  W: 


s 


Oil 

*--1 

12.0000 

46.5882 

67 . 2941“ 

✓ 

^i2 

-1 

9.6613 

23.1534 

0 

^i3 

-1 

6.8313 

0 

0 

^i4 

_  -1 

0 

0 

0 

✓ 

✓ 

~  2k 

and  W  contain  all  of  the  coefficients  of  Q^s  resulting 
from  the  expansion  of  (4-142).  Still  to  be  found  are  the 
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r\  -i  “** 

coefficients  of  not  involving  any  ,  the  o^'s  and  oi*'s. 

Therefore,  from  (4-79)  and  (4-80): 

(s4+12s3+48s2+80s+48) (s4-12s3+48s2-80s+48) 


=  s8  -  48s6  -  480s4  -  1792s2  +  2304, 


a 


a 


a 


0 

1 

2 


2304 

-1792 

480 


and 

(s4+19s3+138s2+435s+450)  x 
(s4-19s3+138s2-435s+450) 


=  s8-85s6  +  3414s4  -  65025s2  +  202500  , 


(4-148) 


(4-149) 


(4-150) 


aQ*  =  202500 
a- L*  =  -65025 
a2*  =  3414 


(4-151) 
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Determination  of  the  Q^j's  results  from  equating 

coefficients  of  like  powers  of  s  from  the  expansion  of  the 

determinants  in  equation  (4-142).  The  Q^'s  are  obtained 

in  a  successive  method  (4-65)  by  starting  with  the  s° 

coefficients,  then  continuing  by  equating  coefficients  of 
2k 

s  ,  ke[l,  2,  ...  n-13,  in  scandent  fashion. 

For  s^: 

from  G,  the  only  non-zero  element  *s  g14>  corresponding  to 
row  1,  therefore,  level  1  in  W.  The  only  non-zero  element  in 
level  1  of  W  in  the  s°  column  is  W-j.41'  Equating  coefficients 
yields : 


<g14  x  W141J  Q11  +  a0  s  a0* 


(4-152) 


Solving: 


202500-2304 


^11  - 

T67T21 

341T1 

[67.2941) 

from 

(4-67) 

and 

(4-68) , 

5l2  = 

*21 

- 

-88.95327 

^13  = 

53i 

= 

48.14819 

^14  = 

Q41 

= 

-3.40294 

44,20803 


(4-153) 


(4-154) 
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1X3  t 


Solving : 


Q33  =  351.24159 


and  from  (4-68), 


Q34  =  Q43  =  -135.68816 


for  coefficients  of  s  : 


4  4- 

E  g4i  Z  W,-ln  Q-i  i  +  =  <*** 


*-  &4X  "il-i  y-i  i 

j=l  31  i=l  113  ]1 


Q44 


4  4 

Z  Z  Q..  +  a-*  -  o 

3*1  1=1  13 

i+ j*8 


*  °3*  -  “3 


<5  =  120.37773 


is  now  entirely  known. 


(4-159) 


(4-160) 


(4-161) 


(4-162) 


(4-163) 


(4-164) 
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44.20803 


-88.95327 


48.14819 


-3.40294 


Q 


-88,95327 

48.14819 

-3.40294 


296.94152 

-263.70162 

55.71337 


-263.70162 
351.24159 
-135 . 68816 


From  equation  (4-91): 


55.71337 

-135.68816 

120.37773 

(4-165) 


Q=  PTQP 


200196 

0 

0 

0 

0 

63233 

0 

0 

0 

0 

2934 

0 

0 

0 

0 

37 

(4-166) 


Substituting  matrices  Q,  R,  A,  B  and  C  into  the  Riccati 
equation,  and  solving,  yields  the  optimal  feedback  control 
law,  u*; 


i*  =  -[402  355  90  7] 


(4-167) 


where  G*  =  -[402  355  90  7], 


(4-168) 


The  optimal  closed  loop  system,  x  =  (A+BG*)x, 
in  phase  variable  form  is: 
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r 


(4-169) 


-450  -435  -138  -19 


which  realizes  the  desired  set  of  eigenvalues. 

3 .  Higher  Order  Example  (Seventh) 

The  previous  examples  represent  an  odd  and  an  even 
ordered  system,  illustrating  the  minor  differences  in  com¬ 
putational  procedures .  It  is  observed  that  for  low  order 
systems,  the  calculations  can  be  done,  relatively  easily  by 
hand.  Higher  order  systems  require,  laborious  and  tedious 
computations.  Appendix  C  provides  a  digital  computer  program 
which  yields  the  weighting  matrix,  Q,  with  the  only  required 
input  being  the  transfer  functions  of  the  known  and  desired 
eigenvalue  systems. 

The  following  seventh  order  example  utilizes  the 
results  from  the  program  given  in  Appendix  C. 


Consider  the  linear  system  given  by  its  transfer 


function : 


T?(s) 


249.435788 


s  +9.0s°  +  4  0.4s  +116. 8s  +233 .6s  +32  3 . 2s ‘ 


+  288.0s  +  128 . 0  . 


(4-170) 
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It  is  known  that  the  characteristics  of  the  desired  system 
are  such  that  its  resulting  transfer  function  is  given  by: 


T?*(s) 


249.435788 

T 


s^+15.4s6+101.64s5+372.68s4+819 .896s3 


(4-171) 


+10 82. 26272S+793. 659328 s+ 249. 4357888 


The  diagonal  Q  matrix  was  determined  from  the  computer  pro¬ 
gram  yielding: 


Oil 

45834.21273428 

Q22 

= 

89780.21841734 

Q33 

= 

55970.39786199 

Q44 

= 

19170.7157376 

^55 

= 

3959.33664 

Q66 

= 

494.9776 

Q77 

= 

33.68 

Q,  R  and  the  state  and  output  matrices  representing  the 
transfer  function  in  equation  (4-170)  in  phase  variable  form 
were  substituted  into  the  matrix  Riccati  equation.  The  opti¬ 
mal  feedback  control  law,  u*,  resulting  from  solution  of  the 
Riccati  equation  is: 
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u*  =  -121.435789  x  -  505.659319  x2 

-759.06269  x3  -  586.29596  x4  -  255.879974  xg 
-61.2399914  xg  -  6.39999944  x?  ,  (4-173) 

where  the  matrix  G*  is: 


G*  =  -[121.435789  505.649319  759.06269  586.29596 

255.879974  61.2399914  6.39999944].  (4-174) 

The  optimal  closed  loop  system, 

x  =  (A  +  BG*)  x 

y  =  Cx  , 

expressed  as  a  transfer  function  is: 


Y(s) 
u*  ( s ) 


249.435788 

T 


- - * - p — 

s  +15. 39999974s  +101. 6399914s  + 


372. 679974s4+ 819. 89596s 3+ 

1082. 26209s2+793.659319s+249. 435789, 


which  realizes  the  desired  transfer  function  with  the  pre¬ 
scribed  eigenvalues. 

Figures  4.2  through  4.13  show  the  impulse  and  step  res¬ 
ponses  of  the  three  previous  examples,  both  before  and 
after  compensation. 
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Figure  4.2.  Unit  Impulse  Response  of  the  Uncompensated  (1) 

and  Prescribed  Eigenvalue  (2)  Third  Order  Systems 


Figure  4.3.  Unit  Step  Response  of  the  Uncompensated  (1)  a 
Prescribed  Eigenvalue  (2)  Third  Order  Systems 
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Figure  4.5.  Unit  Step  Response  of  the  Compensated  (1)  and 
Prescribed  Eigenvalue  (2)  Third  Order  Systems 
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Prescribed  Eigenvalue  (2)  Fourth  Order  Systems 


Figure  4.10.  Unit  Impulse  Response  of  the  Uncompensated  (1) 
Prescribed  Eigenvalue  (2)  Seventh  Order  Systems 


Figure  4.11.  Unit  Step  Response  of  the  Uncompensated  (1)  and 
Prescribed  Eigenvalue  (2)  Seventh  Order  Systems 
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Prescribed  Eigenvalue  (2)  Seventh  Order  Systems 


CONCLUSIONS  AND  DISCUSSION 


A  method  has  been  shown  for  determination  of  the  state 
weighting  matrix  in  order  to  satisfy  a  prescribed  set  of 
eigenvalues  through  phase  variable  state  feedback.  From  a 
strictly  mathematical  viewpoint,  this  technique  requires  only 
a  knowledge  of  matrix  algebra.  Every  attempt  has  been  made 
to  avoid  the  necessity  of  inverting  a  matrix.  The  intro¬ 
duction  of  Chapter  IV  made  known  the  fact  that  previous  devel 
opments  in  this  area  have  suffered  the  main  drawback  of  res¬ 
triction.  The  author  believes  the  method  presented  here, 
using  Cauer  Second  Form,  overcomes  many  of  these  restrictions 
It  presents  a  rational  computational  procedure  for  determina¬ 
tion  of  the  weighting  matrix,  Q;  the  system  eigenvalues  are 
only  required  to  be  in  the  left  half  of  the  complex  plane  as 
opposed  to  the  left  of  a  line  parallel  to  the  imaginary  axis ; 
and  the  method  is  no  more  complicated  for  multiple  or  complex 
eigenvalues  than  a  system  with  linearly  independent  eigen¬ 
values  (or  eigenvectors). 

It  should  be  noted  that  some  authors  "define"  the 
eigenvalue^ )  of  a  matrix  to  be  only  the  real  root(s)  of  the 
characteristic  equation.  In  the  development  that  has  pre¬ 
ceded  in  this  thesis,  all  roots  are  considered  eigenvalues 
of  the  associated  matrix. 


The  algorithm  derived  in  Chapter  IV,  is  designed  as  a 
basis  for  future  research.  In  particular,  if  the  designer 
is  working  with  nth  order  systems  where  n  is  relatively 
large,  it  may  be  advantageous  to  look  at  using  mth  order 
simplified  models  (m<n)  of  each  system  by  a  partitioning 
scheme  similar  to  that  in  Chapter  III.  Again,  it  is  em¬ 
phasized  that  reduced  order  models  do  not  necessarily  yield 
stable  systems.  If  the  simplified  system  retains  the  basic 
characteristics  of  the  original  system,  especially  in  steady 
state,  then  this  would  appear  to  be  a  reasonable  approach. 

A  parallel  approach  could  also  be  inventigated  regarding 
multi-input  multi-output  (MIMO)  systems,  treating  each 
element  of  the  system  transfer  matrix  as  an  individual 
transfer  function. 

To  the  author's  knowledge,  no  work  has  been  done  in  the 
digital  or  sampled  data  areas  involving  continued  fraction 
theory.  This  area  should  be  considered  due  to  the  increasing 
use  and  need  for  digital  control  systems. 

These  topics  represent  just  a  few  of  the  areas  available 
for  future  work. 
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APPENDIX  A 


Consider  the  system 

x  =  Ax  +  Bu,  (A-l) 

where  x  is  an  n-dimensional  state  vector,  u  is  the  input 
function,  and  A  and  B  are  time-invariant  (nxn)  and  (nxl) 
matrices,  respectively.  The  phase  variable  (canonical) 
system  representation  is  defined  as 

A  A 

v  =  Av  +  Bu,  (A- 2) 

where  v  is  an  n-dimensional  state  vector  and 


A  = 


0 

i 

0  .  .  . 

...  0 

0 

0 

0 

1  .... 

...  0 

0 

0 

0 

0 

0 

A 

0 

• 

• 

• 

• 

,  B  = 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

0 

0 

0  .  .  . 

...  1 

• 

« 

■Q.-, 

-a„ 

...  -a 

1 

L 1 

2 

3 

(  A-  3 ) 


The  systems  represented  in  (A-l)  and  (A-2)  are  said  to  be 
equivalent  if  and  only  if  there  exists  a  non-singular  matrix, 
K ,  such  that 


Kv 


(A- 4 ) 
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Kalman  [23]  has  shown  that  a  necessary  and  sufficient  condi¬ 
tion  for  such  an  equivalence  to  exist  is  that  the  system 
in  (A-l)  be  completely  controllable. 

The  controllability  matrix  of  system  (A-l)  is  defined  by 


E  =  [B 


(A-5 ) 


or  in  an  equivalent  manner 

E  =  [e,  e?  ....  e  ],  (A-6) 

~  12  n 

where  the  (nxl)  vector  e^  is  recursively  defined  as 

•i+i  "  6  Si  •  !i  =  5  •  <A-7> 

/V 

The  controllability  matrix  of  system  (A-2),  E,  is  defined 

A  /V 

in  a  similar  manner  with  A  and  B.  Since  there  is  only  one 
control  input,  a  necessary  and  sufficient  condition  for 

A 

controllability  is  that  the  (nxn)  matrix  E  (or  E)  have  an 
inverse . 

Silverman  [20]  has  shown  that  if  the  system  in  (A-l) 
is  controllable,  then  the  transformation  matrix,  K,  is 
determined  by 

K  =  E  E"1  ,  (A-8 ) 


where 


The  elements  of  a  are  the  coefficients  of  the  characteristic 
polynomial : 


det [SI-A]  =  detCSI-A]  =  S 


n  ] 
+  Z  a.S 

i=l  1 


( A-ll ) 


The  matrix  inversion  in  equation  (A-10)  can  be  avoided  by 
using  the  Leverriec-Fadeev  method  for  calculating  the 
coefficients  of  the  characteristic  polynomial.  Once  the 
coefficients  are  known,  E  is  written  by  inspection. 
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Rane  [24]  presented  a  simplified  procedure  for  finding 
the  transformation  matrix,  K,  requiring  no  matrix  inversions. 
Substituting  equation  (A-4)  into  (A-l)  and  premultiplying 
equation  (A-2)  results  in: 

x  =  AKv  +  Bu  ( A-12 ) 

=  KAv  +  KBu  .  (A- 13) 

Comparison  of  equations  (A-12)  and  (A-13)  yields 

AK  =  KA  (A- 14 ) 

and 

B  =  KB  .  (A-15) 

Partition  K  into  n  column  vectors,  each  (nxl),  so  that 

K  =  [k,  k_  _ k  ].  (A-16 ) 

~  ~  1  ~  l  ~n 

Substitution  of  (A-3)  and  (A-16)  into  (A-14)  and  (A-15)  gives 

A[k^  kj  k^  ....  k^ ]  — 

0  1  .  .  .  .  0  0 

0  0  .  .  .  .  0  0 

•  •  •  • 

•  •  •  * 

•  •  •  • 

•  •  •  • 

0  0  .  0  i 

a!  “a2 - an-l  "a 
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B  =  [k.  k_  k~  ....  k  ] 
~  ~  1  ~  2  ~  3  ~n 


From  (A-17 )  and  (A-18) 


k  =  B 
~n 


k  ,  =  A  k  +  k  a 
~n~ 1  ~~n  ~n  n 


k  0  =  A  k  ,  +  k  a  , 
~n- 2  ~  ~n-l  ~  n  n-i 


A  k_  +  k  a. 
-  ~  3  -n  3 


6  *2  *  in  a2  • 


or,  in  general, 


k  .  =  A  k  .  , .  +  k  a_  .  . , 
~n-x  ~  -n-i+1  ~n  n-i  +  1 


k 

-n 


(A-18) 


(A-19 ) 


(A-20 ) 


for  ie[l,  2,  ...,  n-1].  The  column  vectors  k^,  ...,  kn  are 
found  in  a  simple  recursive  manner  and  completely  determine 
the  transformation  matrix. 
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APPENDIX  B 


INVERSION  OF  CAUER  I  AND  CAUER  II  FORMS 

This  program  was  written  in  FORTRAN  IV,  and  requires 
minimal  input.  The  only  information  required  is: 

1.  the  order  of  the  system 

2.  which  inversion  is  required 

3.  the  quotients  from  the  continued  fraction  expansion. 
Multiple  data  sets  are  possible,  and  input  in  the  following 
format : 


Card 

Columns 

Description 

Format 

1 

1-3 

M  =  the  desired  order  transfer 

13 

function 

2 

1-20 

h,,  the  first  quotient  of  either 
Cauer  1  o*  Cauer  II  continued 

D20.13 

fraction  expansion 

21-40 

h2,  the  second  quotient 

D20.13 

41-60 

h3 ,  the  third  quotient 

D20.13 

61-80 

h,, ,  the  fourth  quotient 

D20.13 

• 

• 

• 

• 

• 

N 

1-20 

hUM_7>  the  (4N-7)th  quotient 
such'that  4N-7<2M 

D20.13 

21-40 

h4N.6’  4N-6±2M 

D20.13 

41-60 

h4N-5’  4N-5i2M 

D20 .13 

61-80 

h4N-4>  4N-4<2M 

D20.13 

Four  quotients  per  card  until 

2*M  quotients  have  been  input, 
where  M  is  the  system  order. 

Assume  this  is  the  Lth  card. 

The  (L+l)th  data  card  begins 
the  second  data  set. 
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L+l  1-3  M=system  order  13 

4-6  K=1  for  Cauer  I  and  K=2  for  13 

Cauer  II  inversion 


The  computer  program  has  been  written  to  handle  up  to 
20th  order  systems  (M_<20).  If  it  is  required  to  work  with 
higher  order  systems,  only  one  card  change  must  be  made. 

The  specification  statement  is  modified  to  read: 

REAL *8  A(N) ,  B(N) ,  C(N)/N*0./,  D(N)/N*0./,  DZERO 

where  N  is  an  integer  no  larger  than  999.  This  restriction 
can  be  lifted  by  changing  statement  2  to  read: 

2  FORMATC  2IR) 

where  R  is  the  mantissa  of  log-^CN)  .  The  REAL*8  in  the 
specification  statement  indicates  that  all  following  varia¬ 
bles  and  arrays  are  real  valued  and  in  double  precision. 
Modification  to  either  single  or  extended  precision  would 
require  changes  in  all  format  statements.  If  this  is 
desireable,  the  user  should  consult  references  [25]  and  [26]. 

Execution  time  has  shown  to  be  less  than  .18  seconds 
for  systems  of  order  10  or  less. 


non  noon  non  noon  non 


//  EXEC  FORTCLG 
//FORT  .SYSIN  DO  * 

RE AL*8  A( 20) ,8(20) , C ( 2 01/20*3 ./» D (20) /20*0 ./» DZERO 


8 

C 

8 


...  READ  IN  SYSTEM  0RDER»  &  WHICH  CAUER  INVERSION 

16  READ( 5*1)  M,K 

1  FORMAT 1213) 

...  READ  IN  QUOTIENTS  FROM  C3NT.  FRAC.  EXPANSION  . 

RE ADI  5(2)  (A(  I) ,8(1) ,I*l»M) 

2  FORMAT (  4020 . 1  3) 

...  DETERMINE  IF  INVERSION  IS  CAUER  I  OR  CAUER  II 

IF ( K.EQ. 1 )  GO  T3  10 

...  CAUER  II  INVERSION  ... 

...  INITIALIZATION  ... 

C  ( M  I  =*  B(M) 

Ml  *  M-l 
D(  M 1)  *  1.0 
D(M)  »  A(M)*C(M) 

...  ITERATION  ... 

DO  3  I *1 1  Ml 
L  «  1*1 
K  -  N-L 
DO  4  J«1,L 
KJ  *  K*J 
KL  »  KJ— 1 
KP  a  KJ+l 

C(  K J  I  a  3  (  K  )  *D(  <  J  ) 

IF ( J.NE.L )  C ( KJ )  a  C( KJ ) +C( KP ) 

IF ( KL.EQ. 0)  GO  TO  5 
O(KL)  «  A(K)*C(<L)*0(KJ) 

GO  TO  4 
5  DZERO  =  1.0 
4  CONTINUE 

D(M)  «  A( K+l ) *C( M ) 

3  CONTINUE 
GO  TO  11 

...  CAUER  II  INVERSION  ... 

...  INITIALIZATION  ... 

10  CIM)  «  B(M) 

MM  =  M-l 

D( MM)  «  A (M ) *C( M) 

D(M)  a  1.0 


DO 
IP 
MI 
DO 
MJ 
ML 
MP 
C(MJ) 


ITERATION 

6  I«  1  *  MM 
*  1*1 

»  M-IP 

7  J*l,  IP 
»  MI  *J 

«  MJ-i 
MJ*l 


IF ( ML. EQ. 
D(ML)  «  A 
GO  TO  7 


B  ( M  I )  *0  ( M  J  )  *C  (  M  J ) 
01  GO  TO  8 
(MI  ' 


)  *C(  M J ) *D( ML  ) 
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nonon 


8  DZERO  >1.0 
7  CONTINUE 
6  CONTINUE 


...  WRITE  OUT  TRANSFER  FUNCTION  IN  TERNS  OF 
NUMERATOR  ANO  DENOMINATOR  COEFFICIENTS 
WITH  APPROPRIATE  POWER  OF  S  ... 


11  WR ITE( 6  ,12) 

12  FORMAT (///l IX, • N JMERATOR* ,15X, 'DENOMINATOR*  , 

110X, 'POWER  OF  S'  ) 

WR  IT  E(  6, 1 3)  DZERO, M 

13  FORMAT (///31X,D23.13,10X,  12) 

DO  14  1*1, M 
MM I  *  M-I 

WRITEI6.1 51  C ( I ) , D ( I), MM  I 

14  CONTINUE 

13  FORMAT!// /6X,D20.13,4X,D20. 13 ,10X, 12) 

GO  TO  16 

END 

C 

C  FOR  ACTUAL  RUN  THIS  CARD  IS  /*  IN  COLUMNS  l  AND  2 

C 

//GO.SYSIN  DD  * 

£  DATA  INPJT 


APPENDIX  C 


DETERMINATION  OF  WEIGHTING  MATRIX 
(Q)  FOR  PRESCRIBED  EIGENVALUES 

This  FORTRAN  IV  program  was  used  exclusively  on  the 
Naval  Postgraduate  School's  IBM  360/67  digital  computer 
and  includes  the  associated  job  control  language  statements. 
The  program  consists  of  a  main  program  and  nine  subroutine 
subprograms.  The  purpose  of  each  subroutine  is  delineated 
below. 


SUBROUTINE 


DESCRIPTION 


READ 

read  in  coefficients  of  numerator  and  denomi¬ 
nator  polynomials  of  both  transfer  functions, 
and  places  each  system  in  phase  variable  form. 

RAMAT 

determines  the  Routh  array  matrix  and  the 
transformation  matrix,  P. 

MULT PH 

multiples  two  matrices,  Y  and  Z,  and  gives 
the  resulting  matrix  YZ . ~ 

HMATRX 

determines  the  quotients  of  continued  fraction 
expansion,  H1(H2),  and  the  state  matrix  in 
Cauer  II  form,'’HHl(HH2)  . 

POLYNM 

T 

determines  the  product  det(SI-A)xdet(SI+A. ) 

HELP 

computes  the  matrices  G,  T,  and  W  as  given 
in  Chapter  IV. 

QTILDA 

computes  the  matrices,  Q,  and  Q  from  results 
of  subroutine  HELP. 

QPIJ 

determines  the_off  diagonal  elements,  q^. , 
of  the  matrix  Q.  ^ 

WRITE 

writes  all  two-dimensional  matrices. 
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The  input  required  has  been  reduced  to  a  minimum.  Multiple 
data  sets  are  possible;  and  input  as  indicated  below: 


Card  Columns 


Description 


Format 


1 

1-3 

N;  the  system  order 

13 

2 

1-16, 

17-32, 

33-48 

a  ,  for  ie[l, 2 , — ,n] ; 

the ^denominator  coefficients  of 

the  known  transfer  function. 

F16 . 8 

• 

•  •  • 

for  an  Nth  order  system,  L 
cards  are  required,  where  L 
=K+1  and  K  is  the  integral  part 
of  N/5  . 

L+2 

1-16, 

17-32, 

33-48 

b  for  ie [ 1 , 2 , . . . , n] ;  the 

numerator  coefficients  of  the 
known  transfer  function,  L 
cards  required. 

F16 . 8 

2L+2 

1-16, 
17-32, 
33-48  , 

•  •  • 

a  ,  for  ie[l,2, . . . ,n]  ; 
tne ^denominator  coefficients 
of  the  transfer  function  with 
prescribed  eigenvalues. 

F16.8 

3L+2 

1-16, 
17-32, 
33-48  , 

•  «  • 

B  ,  for  ie[l, 2 , . . . ,n] ;  the 
numerator  coefficients  of  the 
transfer  function  with  pre¬ 
scribed  eigenvalues. 

F16.8 

for  multiple  data  sets,  repeat 
the  same  prodecure.  Each  data 
set  requires  4L+1  cards. 


This  program  has  been  written  to  accept  systems  up  through 
20th  order.  To  increase  the  capability  of  the  program,  only 
the  dimension  statements  and  the  second  continuation  card  of 
the  equivalence  statement  require  modification.  The  system 
order  capability  can  be  increased  to  50.  Beyond  50th  order, 
the  program  requires  an  excessive  amount  of  storage  space 
(>510K  bytes).  Even  this  limitation  is  easily  overcome  by 
removing  the  four  cards  between  statements  1000  and  1001. 


i 
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Any  other  modifications  (i.e.,  single  or  extended  pre¬ 
cision)  require  extensive  changes  to  all  subprograms.  In 
this  case,  the  user  should  consult  [25]  and/or  [26].  It  is 
recommended  that  an  object  deck  or  disk  storage  be  used  when 
available  as  compilation  time  is  approximately  70-80%  of 
total  CPU  time. 
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CALL  READ(N,AA,BB,CC»DE) 

CALL  RAMAT l N» RAM2  »P2  *  CC  »0E»  R 
CALL  HMATRX(N,RAC0L2 , H2.HH2 » 


CALL  WRITE? N*HH2 ) 
CALL  POLYNMlNtOE, E2I 
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*  WWW 


SUBROUTINE  R£A0(N«Al»Bl#CltOl)  ,  „ 

IMPLICIT  REAL+8  <A-H,0-2J,  INTEGER  (I-N,$) 

DIMENSION  A1I20  «20 )  t  BM20),  CK20J,  01(20) 
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TURN 


SUBROUTINE  RAMAT ( N,R, P,CR,DR,RC ) 

IMPLICIT  REAL*8  <A-H»0-Z),  INTEGER  <I-N,$I 

DIMENSION  R{42,21»»  CR(20»,  DR(20J,  P(20,20),  RC(42) 
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SUBROUTINE  HELP(N) 

IMPLICIT  REAL*8  (A-H,0-Zl,  INTEGER  C  I-M  « S  > 

DIMENSION  H(20f  20) •  AI0(20,20)t  V<20,2 3),  F(20,20,20» 
COMMON  H/THREE/  V.  F 
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RETURN 
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SUBROUTINE  WRITE (N»RHO) 

IMPLICIT  RE  AL*3  <A-H,0-Z>»  INTEGER  U-N.JI 
01 MENS I  ON  RH0(20,20I 
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